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ABSTRACT 

We introduce Copernicus Complexio (COCO), a high-resolution cosmological N-body simu¬ 
lation of structure formation in the ACDM model. COCO follows an approximately spherical 
region of radius ~ 17.4/i -1 Mpc embedded in a much larger periodic cube that is followed at 
lower resolution. The high resolution volume has a particle mass of 1.135 x lO 5 ft- 1 M 0 (60 
times higher than the Millennium-II simulation). COCO gives the dark matter halo mass func¬ 
tion over eight orders of magnitude in halo mass; it forms ~ 60 haloes of galactic size, each 
resolved with about 10 million particles. We confirm the power-law character of the subhalo 
mass function, N(> p) ex p~ s , down to a reduced subhalo mass M su j,/M 2 oo = p = 10 -6 , 
with a best-fit power-law index, s = 0.94, for hosts of mass (M 200 ) = 10 12 /i _1 Mq. The 
concentration-mass relation of COCO haloes deviates from a single power law for masses 
M 200 < a few x 10 8 /i _1 Mq, where it flattens, in agreement with results by Sanchez-Conde et 
al. The host mass invariance of the reduced maximum circular velocity function of subhaloes, 
v = V max /V 200 , hinted at in previous simulations, is clearly demonstrated over five orders of 
magnitude in host mass. Similarly, we find that the average, normalised radial distribution of 
subhaloes is approximately universal (i.e. independent of subhalo mass), as previously sug¬ 
gested by the Aquarius simulations of individual haloes. Finally, we find that at fixed physical 
subhalo size, subhaloes in lower mass hosts typically have lower central densities than those 
in higher mass hosts. 
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1 INTRODUCTION 

Since its introduction over thirty years ago, the cold dark matter 
(CDM) model of structure formation d Peehies|| 1 982[ |Davis et al.| 
1 1985| [Bardeen et al.|19 86) has been extensively investigated theo¬ 
retically and tested with an impressive array of observational data. 
According to this, the now standard, model of cosmogony, galaxy 
formation is driven by the evolution of the dark matter (DM) haloes 
in which the galaxies reside. It is into these haloes that gas cools 
and condenses, becomes unstable and fragments into stars, leading 
to the formation of galaxies ( [White & Rees|1978| |White & Frenk| 
|1991| >. This basic picture has been elaborated in detail using simula¬ 
tions and semi-analytic models and it has largely been confirmed by 
countless observations (see e.g. |Frenk & White||2012[ for a recent 
review). Thus, DM haloes are the fundamental non-linear building 
blocks of cosmic structure [Kauffmann, White & Guiderdoni|1993| 
|Cole & Laceyl 199~6[ and understanding their properties, abundance 
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and spatial distribution has been a subject of extensive theoretical 
study over the past thirty years. 


The formation and evolution of DM haloes and of the galaxies 
residing within them is modelled in the context of the background 
cosmological model that describes the expansion and growth his¬ 
tory of the Universe. The last 20 years have seen the emergence 
of the ’’Lambda Cold Dark Matter” (ACDM) model, which com¬ 
bines a flat Friedmann-Lemaitre model with a cosmological con¬ 
stant - A -, responsible for the late-time accelerated expansion of 
the Universe, with the CDM paradigm, in which thermally cold 
relic elementary particles constitute the majority of non-relativistic 
matter and govern the growth of cosmic structures. The predictions 
of the ACDM model have been tested to a high precision on lin¬ 
ear scales, from the very early universe {e.g. [Hinshaw et al.|20T3} 
Planck CollaborationetaL|2014 2015b to large cosmic scales (e.g. 
Cole et al.|2005||Eisenstein et al. 12005]|Percival et al.|2010||Ander-| 


et al,|2012 |de la Torre et al.|20T3 |Driver et al.|2009[ , yielding 


very good agreement with observations. To further test and con¬ 
strain the current model, one needs to study its predictions down to 




























2 Wojciech A. Hellwing et al. 


smaller scales, extending significantly into the non-linear regime of 
structure formation. 

N-body simulations represent the most widely used and con¬ 
venient method of exploring the highly non-linear regime of cos¬ 
mic structure formation. Starting from a set of initial conditions, 
the numerical simulations follow the formation and evolution of 
structures from an early epoch down to present day. Motivated by 
the fact that DM represents most of the matter in the Universe and 
because of the relatively simple physics of collisionless DM par¬ 
ticles, DM-only simulations represent the most widely used cate¬ 
gory of numerical simulations. When designing a cosmological N- 
body experiment, one is concerned by two major factors. Ideally, 
one would like to simulate a region of the universe that is as large 
as possible to get a representative census of the structures encom¬ 
passed within it. On the other hand, one would also want very high 
mass resolution, to be able to resolve accurately even the smallest 
cosmologically relevant objects. Unfortunately, due to limited com¬ 
putational resources, these two requirements are in conflict, which 
implies that various compromises need to be made when designing 
a numerical simulation. So far, the biggest efforts were focused into 
two, somewhat complementary approaches. The first is represented 
by simulations like Millennium (MS) jSpringel et al,||2005) , Mil¬ 
lennium II (MS-II) l Boylan-Kolchin et al.|2009 |, Millennium XXL 
(MXXL) l Angulo et al.|20 1 2|, BolshoTi Klypin, Trujillo-Gomez &| 


|Primack|2011) , MultiDark |Prada et al.|20 1 2\ . Horizon Run I-III 

(|Kim et al.|2009[|2011[ >, Horizon-47r jPrunet et al.|2008 |Teyssier| 
|et al.|2009| >, MareNostrum Universe ( [Gottloeber et al.|2006} , Ju- 
bilee project (| Watson et al.|2014) , Coyote Universe jHeitmann et al.| 
|2010| >, DEUS simulation ( |Alimi et al,||2012| |Rasera et al,||2014) 
or MICE suite (Fosalba et al.||2015} . These follow stmcture for¬ 
mation in a large cosmological volume at the expense of having 
a medium or a low mass resolution. Such simulations provide the 
formation histories for a very large number of medium- and high- 
mass DM haloes, but do not necessary resolve all the details rele¬ 
vant for galaxy formation. On the other side we have N-body simu¬ 
lations like the AQUARIUS project i Springel et al.|[2008} , the Via 
Lactea (jDiemand, Kuhlen & Madau i2007J, the Phoenix project 
(|Gao et al.||20 1 2|, CLUES (|Gottloeber, Hoffman & Yepes||20 1 0) 
and the ELVIS suite l |Garrison-Kimmel et al.|2014b[ l that are char¬ 
acterised by a very high mass and force resolution but are limited 
to very small cosmic volumes. These give a very detailed picture 
of galaxy- and cluster-size haloes, but do so only for a very lim¬ 
ited number of objects, which makes their results sensitive to small 
number statistics, and are unable to capture the full interconnec¬ 
tion between small (DM haloes) and large (the cosmic web) cosmic 
scales. 

The recent years have seen a lot of attention focused on ob¬ 
taining detailed histories for a large number of Galatic-size DM 
haloes. This is because our own Milky Way (MW) Galaxy together 
with the Local Group (LG) galaxies, thanks to their direct prox¬ 
imity, constitute an important test-bed for cosmic structure forma¬ 
tion theories. Thanks to an ever growing accuracy of astronomical 
observations, we are presented with a very detailed picture of our 
nearest cosmic neighbourhood. The past decade has brought an im¬ 
pressive amount of data on the MW, Andromeda and their satellites, 
as well as on other small members of the LG ( e.g. Belokurov et al.| 


2006 2007 2014} McConnachie et al.|2009|[KoposoT^raL|2008[ 


McConnachie||20 1 l\ . These data have led to a number of appar¬ 

ent discrepancies between the predictions of numerical simulations 
and observations, which are collectively known as the “ACDM 


small-scale crisis”. The “missing satellites problem” (Kauffmann. 


|White & Guiderdoni| 1993||Klypin et al.|1999b||Moore et al.| 1 999^ 


was among the first to be recognised. Here the tension arises due 
to the fact that dissipationless numerical simulations predict many 
more small DM satellites (clumps or subhaloes) in a Galactic-size 
halo than the actual number of observed MW satellites. One of the 
most favoured solution to this problem predicts that below a cer¬ 
tain mass-scale the majority of DM satellites have failed to host 
luminous galaxies. This is due to baryonic physics, related to hy¬ 
drodynamic, energy feedback processes and the reionisation, that 
depletes the cold gas from small mass haloes, thus preventing star 
formation and rendering these objects dark (for the most recent re¬ 
sults see e.g. |Schaye et al.||2015| |Boylan-Kolchin||20 1 4| | Vogels- 1 
Iberger et al.|2014||Sawala et al.|2015||2014| >. 

Another small-scale ACDM discrepancy, emphasised in re¬ 
cent years by |Boylan-Kolchin, Bullock & Kaplinghat|(201 1| , is the 
so-called “Too Big Too Fail” problem. It is due to the inconsis¬ 
tency between the internal kinematics of the observed 11 classical 
dwarf MW satellites and the distribution of kinematic parameters 
inferred for the most massive satellites of MW-size hosts in the 
AQUARIUS simulation suite i jBoylan-Kolchin, Bullock & Kapling^j 
|hat||20l2|. Recently, a similar claim was made also for the field 
dwarf galaxies found in the LG ( jGarrison-Kimmel et al.|]2014a[ l. 
This discrepancy has various possible solutions, being a possible 
manifestation of highly non-linear and stochastic baryonic physics 
in low mass haloes jSawala et al. ||20 16f and the impact of stel¬ 
lar feedback on DM density profiles (see e.g. Pontzen & Goveu| 
|nato|2012||Ofiorbe et al.|2015HBrook & Di Cintio|20 1 5| i. Others 

have also shown that the problem can be largely alleviated when 
the mass of the MW (LG) is sufficiently low (e.g. |Wang et al.|2012[ 
ICautun et al.|2014a] >. 

Finally, the recent discovery that a subset of Andromeda satel¬ 
lites are distributed in a thin plane (e.g. |McConnachie & Irwin| 
|2006[|Ibata et al.|2013| and their radial velocity components show 
some degree of a coherent co-rotation, together with the previ¬ 
ously known thin polar disk-like distribution of the MW satellites 
(e.g. |Metz, Kroupa & Libeskind|2008HPawlowski & Kroupa|2013| >, 
were postulated by some authors to also present a challenge for the 
ACDM paradigm. 

It is important to note that many of these apparent points of 
tensions were derived from comparisons with a rather small num¬ 
ber of host haloes. In particular, obtaining sufficient resolution to 
study the internal properties of subhaloes in a MW-size host ne¬ 
cessitates the use of ultra-high resolution simulations, which are 
limited to very small volumes and only a handful of central host 
haloes. Given such a limited sample of host haloes, the resulting 
satellite populations may be prone to halo-to-halo scatter, which is 
intrinsic to hierarchical models for structure formation. MW-size 
haloes are characterised by variety of evolutionary histories and 
large-scale structure environments in which these systems evolve. 
These are important factors that need to be properly evaluated and 
understood before claiming any potential discrepancies between 
ACDM galactic-scale predictions and the MW and LG observa¬ 
tions. For example, proper cosmological-volume simulations could 
help us determine to which extent our own Galaxy, its DM halo and 
the satellite system are rare or special within the ACDM paradigm. 
This concerns one of the core assumptions of modern cosmology, 
the Copernican Principle , according to which a MW-based ob¬ 
server is not privileged in the sense that we can observe a fair sam¬ 
ple of the Universe. 

In this paper we introduce the Copernicus Complexio (the 
Copernicus Conundrum; hereafter COCO), which is a DM-only 
simulation tailored for the study of a statistically significant sam¬ 
ple of well resolved MW-size haloes and their satellites. The 
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Figure 1. Projected density along the z-axisof a 70.4 X 70.4x 1.5/i _1 Mpc 
slice centred on the middle of the COCO simulation at redshift 2 = 0. The 
various colours show the density at different resolution levels: lowest res¬ 
olution (grey), medium resolution (orange and purple) and high-resolution 
(blue to yellow). Note the amazing level of the cosmic web details seen 
inside the high-resolution region. 

simulation follows a hybrid ’zoom-in’ approach, similar to the 
one adopted in the GIMIC simulation suitel jCrain et al.||2009f 
(Galaxies-intergalactic medium interaction calculation), with a 
high-resolution region of radius ~17/F 1 Mpc embedded within 
a much larger box resolved at low-resolution. The large volume of 
the high-resolution region contains around 60 MW-size haloes and 
their satellite populations, resolved at a resolution close to that of 
the AQUARIUS level 3 simulations. This is more than sufficient to 
properly capture the internal structure and properties of subhaloes 
hosting faint MW satellites, attaining at the same time a good statis¬ 
tical sample of DM hosts of various masses located in diverse envi¬ 
ronments. In addition, the simulation contains a very large number 
of well-resolved lower mass haloes, whose properties are studied 
here for the first time with such good statistics. 

In this paper we introduce the new COCO simulation and 
present the first-stage analysis of its results. Section 2 presents our 
selection of the high-resolution region and gives details on the nu¬ 
merical and cosmological set-up. The results on DM halo abun¬ 
dances, formation times and internal density profiles are presented 
in Section 3. In Section 4 we study in detail the populations of 
satellite subhaloes, including their mass and velocity functions, ra¬ 
dial distributions, internal kinematics and effects of host-induced 
tidal stripping. We give our concluding remarks in the Section 5. 


2 THE COPERNICUS COMPLEXIO COSMOLOGICAL 
SIMULATION 

The COCO simulation was designed with the goal of resolving the 
formation and evolution of MW-size haloes and their subhaloes in 
a representative cosmological volume. This prompted the use of a 
zoom-in simulation (Katz & White|1993l|Frenk et al.|1996l|Crain| 


et al. 2009; Ofiorbe et al. 2014) that captures in very great detail the 
evolution of a selected region, which in turn is embedded within 
a larger cosmological volume that is simulated at low-resolution. 
The role of the latter is to produce the correct large scale modes 
and tidal fields inside the high-resolution region. Starting from a 
low-resolution simulation, the high-resolution volume was selected 
by optimizing the number of Galactic-mass haloes that could be 
resolved given the available computational resources. 

Randomly selecting the high-resolution region can result in it 
containing one or more rich clusters. Such massive objects would 
dominate the computational time required for the whole COCO sim¬ 
ulation, leading to a wastage of resources, since we are primarily 
interested in MW and lower mass objects. To avoid unnecessary 
computations, but keeping in mind that we want to simulate a fair- 
sample of the Universe, possibly close to the observed Local Vol¬ 
ume, we have selected a region that satisfies the following criteria: 

(i) there are no cluster-mass haloes (M it 5 x 10 13 /i _ 1 Mq) 
inside the zoom-in region, 

(ii) there are no massive cluster haloes (M > 5 x 10 14 /i _ 1 Mq) 
within 5/z 1 Mpc of the zoom-in boundary, 

(iii) the mass function of MW-mass haloes (M ~ 10 12 /i _ 1 Mq) 
is as close as possible to the universal mass function. 

2.1 Cosmological and numerical parameters 

The COCO simulation follows structure formation in a high- 
resolution region that is approximately a sphere of radius ~ 
17.4/i _1 Mpc (2.2 x 10 4 ft -3 Mpc 3 in volume) embedded within a 
70.4/t -1 Mpc low-resolution periodic box, as illustrated in Fig. |T] 
It uses a Wilkinson Microwave Anisotropy Probe (WMAP) - sev¬ 
enth year result cosmogony ( |Komatsu et al.|20TT) with the follow¬ 
ing cosmological parameters: 

D m0 = 0.272 , Gao = 0.728 , Q b = 0.04455 , 

= 0 , h = 0.704 , cr 8 = 0.81, n s = 0.967. (1) 

The high resolution region consists of ~12.9 billion particles of 
mass 1.135 x 10 5 /i _ 1 Mq and of ~510 million medium and 
low resolution particles that have progressively larger masses. The 
Plummer equivalent force softening was chosen to increase from a 
value of 0.23/i -1 kpc for the high resolution particles to a value of 
23/i _1 kpc for the lowest resolution level. 

The high-resolution region was selected from a lower resolu¬ 
tion version of the COCO simulations that we refer to as COpernicus 
complexio LOw Resolution (COLOR). The COLOR simulation has 
the same corresponding initial phases as COCO but is set-up with 
1620 3 DM particles uniformly distributed throughout the whole 
70.4/t -1 Mpc periodic box. It has a mass and force resolution of 
~ 6.2 x 1 O 6 /i - 1 M 0 and l/i _1 kpc respectively, which is exactly 
the same as in the MS-II. While the COLOR box size is relatively 
small, the lack of very large-scale modes has little impact on the 
internal properties of galactic and smaller mass haloes (for more 
details see jPower & Knebe|2006| >. 

The selection of the COCO zoom-in region was performed by 
generating a large number of randomly placed spheres of radius 
17/i -1 Mpc in the COLOR volume at z = 0. We discarded any 
volume not fulfilling criteria (i) and (ii) given in Section [2] From 
the remaining volumes, we selected the one whose halo mass func¬ 
tion in the range, M < 1 O 12 /! _ 1 M 0 , showed the closest match to 
the universal mass function. Following this, the initial conditions 
of the zoom-in region were generated using the same method as 
in the AQUARIUS and the GIMIC JCrain et al.|2009| projects. The 
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particles from the selected z = 0 volume were traced back to their 
Lagrangian positions. The Lagrangian volume was divided in 256 3 
regular cells and each cell occupied by one or more of these par¬ 
ticles was classified as high resolution. The remaining cells were 
classified as medium or low resolution cells depending on the dis¬ 
tance to the nearest high resolution cell. Each high resolution cell 
is filled with a periodic glass distribution of 24 3 particles, while the 
medium to low resolution cells were sampled with progressively 
fewer particles. Higher frequency power was added to the resulting 
particles down to the Nyquist frequency while making sure that the 
lower frequency modes were the same as in the COLOR simulation. 

The initial conditions (initial positions and velocities of all 
particles) for the COCO simulation were set at 2 = 127 using 
second order Lagrangian perturbation theory using the method of 
Jenkins (2010). The initial phases for both the COCO and COLOR 
are taken from the public multi-scale Gaussian white noise field 
called Panphasia, and are published in Table 6 of |Jenkins| \2() 1 3} 
under the alternative name of the ‘DOVE’ simulation. The COCO 
initial conditions differ from COLOR by a uniform spatial trans¬ 
lation so that the coordinate origin in COCO is located at the co¬ 
ordinates (7,16, 44)/t _1 Mpc within the COLOR simulation. This 
translation places the high resolution region of COCO at the centre 
of the simulation volume. 

The COCO and COLOR simulations were both run with GAD¬ 
GETS Tree-PM N-body code, which is an updated version of the 
publicly available GADGET2 code ( |Springel|[2005| l. GADGET3 is 
a hybrid code in which the long-range forces are computed using 
a particle-mesh method, while short-range forces are obtained by 
using a hierarchical oct-tree algorithm. This particular heteroge¬ 
neous architecture allows for a relatively easy follow-up of nested 
grids placed with increasing accuracy around the high-resolution 
region. This results in a proper long-range force accuracy through¬ 
out the box, while focusing most of the computational effort in¬ 
side the high-resolution region of interest. For both simulations DM 
particle positions and velocities were saved in 160 equally spaced 
in log(z + 1) snapshots. Table [I] summarizes some details of the 
COCO and COLOR simulations and compares these with some other 
widely used cosmological simulations. 


2.2 Halo and subhalo finding 

We identified DM haloes and subhaloes using the SUBFIND algo¬ 
rithm (Springel et al .|2001| ). Due to the large number of particles 
and high clustering level of the COCO simulation, the standard ver¬ 
sion of SUBFIND would have required a vast amount of computer 
memory and CPU time. To overcome this problem, we have used 
an updated version of the algorithm that has been especially op¬ 
timised for parallel computing and big data. While these changes 
significantly decrease the required computational resources, they 
do not affect the final output of the method, with the new SUB- 
FIND version producing the same halo and subhalo catalogues as 
the older version. 

SUBFIND starts by identifying DM haloes using the friends- 
of-friends (FOF) algorithm jPavis et al.|1985[ >, for which we used 
a linking length b = 0.2 times the mean inter particle separation. If 
any resulting FOF group has one or more low resolution particles, 
we exclude it from further analysis since the internal properties of 
such objects might have been affected by unrealistic two-body scat¬ 
tering and self gravity. All pristine FOF groups with at least 20 
particles were kept for further analysis. At 2 = 0 we found more 
than 1.2 x 10 7 (5.18 x 10 6 ) FOF groups in the COCO (COLOR) 
run, with the peak value was found at 2 = 3 and consisted of 


1.63 x 10 7 (6.19 x 10 6 ) groups. SUBFIND further analyses each 
FOF group to find gravitationally self-bound DM subhaloes (i.e. 
substructures within the FOF groups). Subhalo candidates are first 
identified by looking for overdense regions inside the FOF groups 
that are further pruned by checking which ones are gravitationally 
self-bounded objects. This results in a catalogue of self-bounded 
structures containing at least 20 particles. For each subhalo we also 
compute and store a number of additional properties. This consists 
of peak circular velocity, V ma x, and the physical radius, R m ax, 
at which this peak is attained, half-mass radius, spin (angular mo¬ 
mentum), position (corresponding to the minimum of gravitational 
potential) and bulk velocity of the subhalo. 

Each subhalo is characterised by a well defined mass, M su b, 
and radius, r au b- The former is given by the mass contained in all 
the particles that pertain to the subhalo. The latter is approximately 
the subhalo proper tidal radius (see Figure 15 of |Springel et al.| 
( |2008} ). The FOF groups are characterised in terms of their FOF 
mass, Mfof, as well as of their M 200 mass. The first, similarly 
to subhaloes, is given by the mass contained in all the particles 
associated to a given FOF group. In contrast, M 200 is the mass 
contained in a sphere of radius r 200 centred on the FOF group, 
such that the average overdensity inside the sphere is 200 times 
the critical closure density, p c . We refer the reader to |Sawala et al.| 
< |2013) > for a comparison of systematic differences between the two 
as well as other halo mass definitions. 

To compute the radial profiles of haloes, we follow a prescrip¬ 
tion similar to that employed by |Power et al.| l |2003) and Knollmann 
& Knebe (20091. Namely, we identify the centre of mass of FOF 
groups using an iterative procedure, by computing the centre of 
mass inside smaller and smaller spheres, with each such sphere cen¬ 
tred on the centre of mass found in the previous iteration step. The 
centre of each FOF group is used to grow logarithmically spaced 
spherical shell bins up to r 200 - Figure[2]illustrates the level of de¬ 
tail to which we resolve MW-mass FOF haloes. 


2.3 Merger Trees 

In CDM cosmologies the first objects to from are DM clumps 
(haloes) with Jeans mass of the order of Earth mass ~ 10~ 6 Mq 
(see e.g. |Green, Hofmann & Schwarz|2004) . Due to numerical lim¬ 
itations, such small density perturbations are not resolved in our 
simulations, hence the first objects to from in COCO have masses 
of ~ 10 e /i -1 MQ, some 12 orders of magnitude larger. However, it 
is wel l established (e.g. |Lacey & Cole|1993[|Kauffmann & White] 
|1993| |Roukema et al.|1997| l that in hierarchical cosmologies char¬ 

acterised by a nearly scale-free Harrison-Zeldovitch like (|Harrison| 
|1970[ |Zeldovich|| 197 2| > initial power spectrum, such as ACDM, 
larger objects forms by consecutive merging of smaller ones. Suc¬ 
cessive populations of haloes grow from mergers of earlier pop¬ 
ulations accompanied by accretion of some smooth mass compo¬ 
nent (see e.g. |Wechsler et al.|2002j l. In order to trace the temporal 
evolution of haloes we constructed DM haloes merger trees (for 
more details see |Helly et al.||2003| . For this, we employed a re¬ 
cently updated algorithm that has been developed for use with the 
semi-analytic galaxy formation code GALFOR M (|Cole et al.|20Q0| . 
The method we used is described in detail inJJimtg^et^L (2014) 
and is an upgrade over the earlier version of |Merson et al. ( |2013| l. 
The essential part of the algorithm consists of unique linking be¬ 
tween subhaloes from two consecutive snapshots. This allows for a 
construction of very precise merger trees at the subhalo level. We 
have applied this algorithm to the COCO simulation, resulting in ap- 
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Table 1. Details of the simulations used and described in this work along with a list of some other large simulations frequently used in the literature. Vt> 0 x 
gives the total comoving volume of the simulation, N p is the number of N-body particles, e denotes the Plummer-equivalent force softening, expressed in 
comoving units, m p is the mass resolution. 


Name 

Vbox [h 3 Mpc 3 ] 

N p 

e [h 1 kpc] 

m p [h 1 M 0 ] 

cosmology 



reference 


COCC0 

COLOR 

~2.2 X 10 4 

3.5 x 10 s 

~2344 : P1 

1620 3 

0.23 

1.0 

1.135 x 10 5 
6.19 x 10 6 

WMAP7 

WMAP7 



this work 

this work 


Millennium-II 

1.0 x 10 6 

2160 3 

1.0 

6.89 x 10 6 

WMAP1 


Boy 

lan-Kolchin e 

al. (2 

009 

Millennium 

1.3 x 10 8 

2160 3 

5.0 

8.61 x 10 s 

WMAP1 



Springel et al. 

2005 


AQUARIUS lvl. 

~1.1 x 10 2 

~530 J pl 

0.12 

~5 x 10 4 

WMAP1 



Springel et al. 

2008 


Via Lactea (LRpl 

~5.5 x 10 2 

~402^ 

0.378 

4.11 x 10 5 

WMAP3 


Diemai 

id, Kuhlen & 

Vladau 

2007 

Horizon-47r 

8.0 x 10 9 

4096 3 

7.6 

7.7 x 10 9 

WMAP3 



leyssier et al. 

(20091 

Bolshoi 

1.6 x 10 7 

2048 3 

1.0 

1.35 x 10 s 

combinational 

Klypin, Tr 

ljillo-Gomez 

& Primack 

2011 


Horizon Run 3 

1.3 x 10 12 

7210 3 

150 

2.44 x 10 11 

WMAP5 



Kim etal. (20111 


Jubilee 

2.2 x 10 11 

6000 3 

no data 

7.49 x 10 10 

WMAP5 



Watson et al. 

2014 


MICE 

2.9 x 10 10 

4096 3 

50 

2.93 x 10 10 

WMAP5 



Fosalba et al. 

2015 



a Here we only consider the high-resolution region. 
^ Actual particle number, N p = 12, 876, 807,168. 
c Actual particle number, N p = 148, 285, 000. 
d X-ray clusters+WMAP5 +SN+BAO. 
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Figure 2. A selection of non-linear structures from the COCO simulation. The left- and right-hand panels show the projected DM density for four MW-mass 
FOF haloes. The central panel gives the projected DM density (in a 1.5h ~ 1 Mpc thick slice) in a region that clearly illustrates the rich hierarchy of non-linear 
structures resolved in COCO, from massive haloes (upper-right comer) to cosmic web filaments and voids. 


proximately 1.312 x 10 9 unique subhaloes contained in the merger 
trees. 


3 DARK MATTER HALOES 

In this section we focus on a few key aspects of DM haloes: 
their abundance as a function of mass, their internal structure and 


their formation histories. Understanding the basic properties of DM 
haloes is a key ingredient of any successful galaxy formation the¬ 
ory, since galaxies are formed and evolve inside their host haloes. 
Furthermore, understanding the link between the properties of DM 
haloes and the luminous galaxies that reside within them is cru¬ 
cial for designing and conducting astrophysical tests of the ACDM 
paradigm. 
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Figure 3. Upper panel: FOF mass functions for the COCO and the COLOR 
simulations, on which we superimpose the Sheth-Tormen prediction (solid 
lines). The COCO results and the ST predictions are plotted for a wide range 
of redshifts, from 2 = 9 (red diamonds) to 2 = 0 (purple circles). For 
comparison we also plot the 2 = 0 results from the COLOR run (open 
hexagons) and the Reed et al. prediction (dashed line). The vertical bars in¬ 
dicate Poisson errors. Bottom panel: The COCO and COLOR mass functions 
normalized by the ST prediction at that redshift. The horizontal dashed lines 
mark a 10% difference level. The vertical dashed (dotted) line illustrate the 
COCO (COLOR) halo mass resolution limit. 


3.1 Mass function 


Accurate theoretical predictions for halo mass functions are needed 
for a number of reasons. For example, they are a primary in¬ 
put for modelling galaxy formation, whether it be physically mo¬ 
tivated semi-analytical models (e.g. |Cole et aT]| 1 994| |2000) or 
statistical-based approaches like abundance matching (e.g. Yang,| 
|Mo & van den Boschfl2003j |Guo et al. 12010} . The abundance of 
haloes across cosmic epochs was studied since the early work of 
|Press & Schechter|j 1 974}, with predictions resulting from the ex 


tended excursion set models based on ellipsoidal collapse (Sheth 
&%rormep|2002j or motivated by N-body simulations (e.g. 


Jenk- 


; et al.|2001[|Warren et al.|2006HReed et al.|2007fr . Such models 


were thoughtfully tested in computer simulations, but only in a lim¬ 
ited mass range Mfof > a few x lO s /t _ 1 M 0 . Due to our unique 
COCO and COLOR simulations, we are able to investigate the abun¬ 
dance of FOF haloes down to lower masses and over a wide range, 
spanning 8 decades in halo mass. 

In Figure [3] we compare the present day COCO and COLOR 
halo mass function with the Sheth&Tormen (ST) prediction jSheth 
& Tormen 2002) and with the improvement suggested by 


Reed 


et al. 1 2007 hereafter R07), which was tuned using results of N- 


body simulations. The R07 prediction includes the dependence 
of the halo mass function on the effective power spectral slope 
(n e f f) at the scale of the halo radius. We find good agreement be¬ 


tween the present day COCO and COLOR mass functions all the 
way from the resolution limit of the COLOR simulation, Mfof ~ 
6 x lO s h^ 1 M 0 , up to the most massive objects found in the COCO 
volume, Mfof ~ 1 O 13 /i _ 1 M 0 . As the resolution limit of FOF 
haloes we adopt a minimum threshold of 100 particles. Note that 
this is different from the resolution limit of converged internal 
(sub)halo properties (e.g. Vmax) that we derive in Appendix [A] 
This assures as that the specific choice of the COCO region (see 
f| 2 ]l did not introduce any significant halo abundance bias (scarcity 
or excess) for the range of halo masses that we are interested in. 
We also see a good agreement with the ST and R07 models at 
2 = 0, though both COCO and COLOR predict slightly more low 
mass haloes. This discrepancy is rather small, for example at a halo 
mass of 10 ' /i _1 Mq, where the difference is the largest, ST predicts 
a ~ 11% lower halo abundance, while the R07 result is ~ 18% 
lower. These differences are unlikely to be caused by numerical ef¬ 
fects. First, while it is well known that the FOF algorithm tends to 
overpredict the abundance of poorly resolved haloes, this effect is 
significant only for objects with fewer than 100 particles (Warren 
|et al.|2006| >. 

Figure[3]also shows the time evolution of the COCO halo mass 
function from redshift 2 = 9 till the present day. The result beau¬ 
tifully reflects the well known hierarchical character of DM halo 
build-up, with smaller haloes forming first, which in turn merge 
into bigger and bigger objects. For comparison for each COCO red¬ 
shift data we also plot the ST prediction line. It is clear that the ST 
prognosis fails to match the COCO data for 0.5 ^ 2 $5 5 as it signif¬ 
icantly overpredicts the abundance of objects. Similar results were 
also found by |Klypin, Trujillo-Gomez & Primack[ l |2011| l. Interest¬ 
ingly this discrepancy is largest for the intermediate redshifts of 
1 ;$ 2 < 2, while at 2 = 9 the ST forecast is again in a good agree¬ 
ment with our data. The epoch at which we observe the biggest dis¬ 
crepancy between the COCO data and the ST predictions also hap¬ 
pens to be the epoch at which the halo merger rates are the highest 
(see e.g. |Flopkins et al.|2010| >, hence the mass function of collapsed 
objects experience the most dynamical evolution, which in turn is 
reflected in the failure of the ST forecast. 

To better compare with analytical models for the abundance of 
haloes, it is more convenient to express the halo mass, M, in terms 
of the variable In z)], where cr(M, z) gives the peak mass 

variance at scale M and redshift 2 . This quantity is defined as 

<t 2 (M, 2) = J p( k , z)W 2 (k, M)k 2 Ak , ( 2 ) 

where P(k, 2 ) is the power spectrum of linear density fluctuations 
extrapolated to redshift 2 and W (k, M) is the Fourier transform of 
the top-hat window corresponding to the radius enclosing mass M 
at the mean density of the universe. Now, the halo mass function 
can be written as 


M dn (M * z) 

AM 


P(z) 


din a 1 
dM 


/( CT ). 


(3) 


where p(z) is the mean mass density of the universe at redshift 2 
and f(a) denotes the halo multiplicity function. This latter quan¬ 
tity, /(<r), takes a universal form that is independent of redshift 
(for more details, see| Jenkins e t al.| 2001 [| Reed et al.|2007 [|Ti nker| 
et al. 2008; Angulo et al. 20f2]». In Figure[4]we plot the multiplicity 
function of FOF haloes as a function of In 2 )] at various 

redshifts. Independent of redshift, the data points follow, to a good 
approximation, the universal shape as predicted by both ST and 
Reed et al. formulas. 
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Figure 4. FOF halo multiplicity function, /(<r), as a function of peak mass 
variance, c(M, z), for a range of redshifts. The vertical bars indicate Pois¬ 
son errors. 


3.2 Density profiles and the mass-concentration relation 


Spherically averaged radial density profiles are one of the sim¬ 
plest yet robust characterisations of the internal structure of DM 
haloes. It is well established that for hierarchical cosmologies like 
CDM, the radial density profile of relaxed DM haloes is to a good 
approximation self-similar and can be mapped by a simple bro¬ 
ken power-law formula, the NFW profile ( Navarro, Frenk & White 
[1996] [1997) : 

E^L =__. ( 4 ) 

pcrit (r/r„)(l + r/r s ) 2 


It characterises the profile of any halo by two parameters: a scale 
radius, r s , and a characteristic overdensity, 5 C . Instead of working 
with the r a and S c parameters, it is customary to define the halo 
concentration, C200, as: 


C200 — 


t"200 

5 

r s 


(5) 


with V 2 oo the virial radius of the halo defined in §2.2| Using this 
parametrisation, the NFW profile effectively becomes a one param¬ 
eter fit, since the characteristic overdensity can be expressed as: 


200 C200 

3 In (1 + C200) — C2oo/(l + C200) 


( 6 ) 


The density profile of a halo is also probed by the shape of the 
circular velocity curve, which is: 


Vc{ r ) = ^MEll : (7) 

where M(< r) is the mass contained inside a sphere of radius r 
centred at the halo centre. For a perfectly spherical halo, the circu¬ 
lar velocity, 14(7), is exactly equal to the circular orbital velocity at 
distance r. For well resolved and relaxed haloes, the circular veloc¬ 
ity takes only one maximum value, Vmax, that is attained at radial 
distance, Rmax . Similar to the virial mass, we can define the virial 


circular velocity V 200 = {GMnoo/i' 200 ) 1 ^ 2 ■ The circular velocity, 
V c (r), bears effectively the same information as the halo density 
profile, p(r), but is much less prone to noise because of the integral 
nature of the former. 

The NFW scale radius, r s , gives the radial position at which 
the r 2 p(r) curve attains its maximum, which sometimes is also de¬ 
noted by r_2 = r s . For the majority of DM haloes, the peak of 
the r 2 p(r) curve is relatively broad. This means that, for haloes 
resolved with a relatively small number of particles, the exact lo¬ 
cation r_2 of the peak is uncertain due to the presence of noise. 
This is reflected in the susceptibility of the NFW fit to the radial 
range used for fitting the profile of haloes resolved with fewer than 
a few thousand particles (Navarro et al.||2004| |Prada et aL|2006[ 


|Gao et al.|2008[|Ludlow et al.|2010 l. This is especially prominent 

when profiles of many similar mass haloes are stacked to remove 
halo-to-halo variation due to the presence of substructures. This be¬ 
haviour indicates that fitting NFW profiles to haloes resolved with 
a relatively small number of particles is biased and gives rise to 
an artificial correlation between concentration and halo mass. As 
a solution, |Navarro et al.] (2004) proposed the use of a more flexi¬ 
ble parametrization, that would account for the differences between 
NFW and stacked universal halo profiles. The improved three- 
parameter fitting formula takes a form, in which the logarithmic 
slope of density assumes a single power law: 

dl°g P(r) _ 2 ( r 

dlogr \ r -‘2 

This induces radial density profile of the form: 


(8) 


In (p(r)/p- 2 ) = -(2 /a)[(r/r- 2 ) a - 1] ; 


(9) 


where p -2 is the density at r_ 2 . The additional parameter a is 
called the shape parameter and, for CDM haloes, it typically takes 
values in the range 0.1 — 0.3. This power-law density profile was 
first introduced by |Einasto[ l |l965) to model the density distribution 
of the stellar halo of our own Galaxy. To distinguish this density fit¬ 
ting function from the NFW profile we will refer to it as the Einasto 
profile. The Einasto profile can be characterised in terms of a con¬ 
centration parameter that is given by eqn. § with r s replaced by 
r_ 2 - Both the Einasto r _2 parameter as well as the NFW scale ra¬ 
dius, r„, correspond to the scale at which the logarithmic slope of 
the density profile attains the ’isothermal' value of -2. In addition, 
for a ~ 0.2 the Einasto profile approximates fairly well the NFW 
profile in the fit range. 

In this work we are interested in a statistical description of 
DM haloes concentrations, with emphasis on the relation between 
concentration and halo mass, its variance and its redshift evolu¬ 
tion. To obtain robust measurements, we fit both the NFW and 
Einasto profiles to all the haloes with at least N™ m = 5000 
particles, which for COCO corresponds to a minimum halo mass, 
M%ob n = 5.7 x 10 8 /i -1 Mq. We discuss further down why we 
picked this particular limiting value. The fitting procedure finds the 
parameter values that minimize the merit function 
N birla 

cr 2 (H) = i ln Pi - ln P/«(“)] 2 . (!0) 

2 = 1 


where the vector of fit parameters E = ( r s ,8 c ) and (r_2, P-2, a) 
for the NFW and Einasto fits, respectively. The number of radial 
bins, Nbins, is equally spaced in log(r) and is selected adaptively 
depending on the number of particles, N p , contained in the halo as: 

Nbins = 2ceiling(6.2 log N p — 3.5). (11) 

This gives Nbins = 92 for our most massive halo and Nbins = 40 
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bins for a halo with 5000 particles. Since the bins are equally 
spaced in log(r), the inner bins contain significantly fewer particles 
than the mid-range and outer bins. Hence, as we move towards the 
halo centre, the radial bins become more affected by sampling noise 
and two body scattering effects. This has been studied by Power 
leFaTI ( |2003t , who found that there is a minimum radius below 
which one cannot trust the radial density profile of haloes extracted 
from N-body simulations. This convergence radius is given by the 
inner most bin that fulfils the Power et al. (2003} criterion (see 
eqn. (20) in their paper). We exclude from the fitting all radial bins 
that are below the convergence radius for a given halo. After apply¬ 
ing this convergence criterion, a halo resolved with N= 5000 
particles is left on average with 20 radial bins. Thus, the minimum 
halo mass for which we perform a fit is given by the mass for which 
more than a half of the radial bins pass the convergence test. 

The individual density profile of haloes resolved with fewer 
than N = 5000 particles is very sensitive to the intrinsic 
numerical noise. However, there is still plentiful of information 
that can be extracted from haloes with fewer than N™ 1 " parti¬ 
cles. By stacking many such haloes, the noise of the density pro¬ 
file is significantly reduced. Finding the position of the maximum 
of the r 2 p(r) curve for a such stacked profile gives the r s and 
r_ 2 parameters for the NFW and Einasto profiles, respectively. 
We apply this method to get median concentration of stack pro¬ 
files for haloes with 800 ^ N p ^ 5000. Using this approach, 
we estimate halo concentrations in COCO down to a halo mass of 
A /200 — 9 x 10 7 h 1 Mq. 

The NFW and Einasto profiles represent a good description of 
the radial density profiles for virialised DM haloes, which are in 
equilibrium. Haloes that experienced recent mergers or close en¬ 
counters can be far away from the state of equilibrium and, thus, 
their density profiles are usually not very well described by neither 
NFW nor Einasto profiles. As a consequence, the concentration pa¬ 
rameter derived from fitting unrelaxed haloes is ill defined and at 
best biased low (see e.g. |Neto et al.|2007j|Gao et al.|2008||Ludlow| 
|et al.|2010[ >. To overcome this problem, we remove non-virialised 
haloes, i.e. objects that do not satisfy the following three criteria 
l |Neto et al.|2007) : (i) the fraction of halo mass contained in its re¬ 
solved substructure is f S ub < 0.1, (ii) the displacement between 
the centre of mass and the minimum of the gravitational potential 
cannot exceed 7% of halo’s virial radius, rooo, and (iii) we require 
that the adjusted virial ratio, v = (2 T — E a )/\U\, is v < 1.35. 
Here T and U are the halo’s total kinetic and potential energy, re¬ 
spectively. To account for the fact that real haloes are not isolated 
objects, we include the Chandrasekhar pressure term, E s , which 
quantifies the degree to which a given halo interacts with its sur¬ 
roundings. See |Shaw et al.||2006| > and their eqn. (6) for the defi¬ 
nition and method used to estimate the pressure term, and also see 
|Power, Rnebe & Knollmann| ( j2012[ l for a more detailed discussion 
about the virial ratio of haloes. 

The left-hand panel of Figure [5] shows the median halo con¬ 
centration as a function of halo mass for our two simulations, 
COCO and COLOR. We find a very good agreement between COCO 
and COLOR results for all halo masses up to M 200 ~ 3.5 x 
10 12 /i _1 Mq. Above this mass, due to scarcity of the massive 
haloes, the COCO results are dominated by halo-to-halo scatter¬ 
ing. This is clearly seen from the increasing size of the error bars 
that show the bootstrap errors of the median. The good agreement 
between the two simulations suggests that we can supplement the 
COCO data at the high-mass end by adding the objects from the 
COLOR simulation. We construct such a joint sample and use it for 
obtaining our best-fit for the median C 200 — A /200 relation as dis¬ 


cussed later. Comparing the results for the NFW and Einasto fits, 
we find clear differences between the two. Below a halo mass of 
M 200 < 1 O 11 /i _ 1 M 0 , the difference takes the form of a system¬ 
atic shift, with the slope of the C200 — A/200 relation being similar 
for both profile fits. 

In the right-hand panel of Figure [5] we also show various fits 
to the C200 (A/200) relation with the goal of comparing the accu¬ 
racy of these literature fits with the results of N-body simulations. 
From the set of single power-law fits, the one that best matches 
the data is the fit proposed by [Jiang et al.| (2014} l. This is in very 
good agreement (better than 4%) with both COCO and COLOR data 
down to a halo mass of ~ 2 x 10 10 /i _ 1 Mg. It seems that below 
that mass our data indicate slight flattening of this relation, hence 
change of the slope. Also the fit of [Klypin, Trujillo-Gomez & Pri-| 
mack| { 20 lT} is in a reasonably good agreement with our data, ex¬ 
cept for the most massive objects (~ 5 x 1 O 13 /i _ 1 M 0 ) where it 
predicts a median concentration that is higher from the value of 
84% of our haloes at that mass. The Klypin et al. results were found 
for halo masses and boundaries defined by a spherically averaged 
overdensity A v i r = 360 x f2 m ,o x p c . This definition roughly corre¬ 
sponds to A 100 in our nomenclature, hence, to allow for a compar¬ 
ison of their fit with our data, we have rescaled their c v i r — M v i r 
relation to appropriate equivalent of C200 — A/200 • However, this 
procedure ideally should be performed for each halo separately at 
the particle level, so our rescaling here can only be treated as an 
approximation. The performance of the |Neto et ak| ( [2007} fit is 
also reasonably good down to A/200 ~ lO 12 /i _i M 0 . The fit of 


|Gao et ah] ( |2008| > agrees with our data only for the most massive 
bins and it clearly predicts a different slope of the concentration- 
mass relation. This unavoidably leads to a significant overestima¬ 
tion of the median halo concentration by their fit for haloes be¬ 
low a mass of ~ 2 x 1 O 12 /i _ 1 M 0 . There are two possible sources 
driving this discrepancy. First, the Gao et al. fit is based on the 
MS which uses significantly different values of the cosmological 
parameters. This difference is most notable for the as parameter 
that is 10% higher in the MS than in our simulations. The dif¬ 
ferent cosmology is probably the main reason for the observed 
differences since variations in both Q. m and as have a large im¬ 
pact on halo concentrations I Ludlow et al.|2013[|Dutton & Macci 6 | 


|2014| |Diemer & Kravtsov|2015[ . Secondly, the Gao et al. relation 

is obtained by fitting a relatively limited halo mass range given by 
5 x 10 12 ^ A/ 2 oo/(A - 1 M 0 ) ^ 10 15 . Since the concentration- 
mass relation is not a simple power law, fitting a single power law 
provides a relation that holds only for that mass range 1 Ludlow et al.| 


|2014HSanchez-Conde & Prada|2014||Prada et al.|2012]T 

To emphasise this last point, we also checked the performance 


of the multi power-law median C200(A/200) model of Sanchez 


|Conde &~Pr ada (2014 hereafter SC14), based on a functional form 
proposed by Lavalle et al.|[2008}: 


5 

(C 200 (A/ 200 , z = 0) ) med = Ci X [In (A/ 20 O /h 1 M 0)] 1 . 

i =0 

( 12 ) 

This multi-component fit is claimed to be a much better match for 
the median halo concentration-mass relation due to flattening of 
this relation at small halo masses. Not surprisingly, the SC14 fit 
shows a very good agreement with the N-body results to an accu¬ 
racy better than 10%. However, the SC14 fit systematically over¬ 
predicts the median concentration for our haloes with A /200 < 
1 O 11 /i _ 1 M 0 . Especially for the first four to five least massive bins, 
we can notice the difference of varying slope of the concentration- 
mass relation between the SC 14 fit and our COCO+COLOR sample. 
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M 20 o [M 0 /h] M 200 [M 0 /h] 

Figure 5. Left-hand panel: Concentration-mass relation for relaxed haloes. The open symbols give the median concentration obtained from fitting NFW and 
Einasto profiles. The best-fitting line was obtained from the combined Einasto COCO+COLOR data, with the COCO concentrations (circles) extended by the 
COLOR results (triangles) for halo masses M200 > 5 X 10 11 h~ 1 M.Q. The filled circles show the average concentrations obtained from stacking profiles for 
the three lowest mass bins. For comparison we also show the data obtain by fitting the NFW profile (open boxes). Right-hand panel: Comparison of various 
C 200 — M 200 fits from the literature (lines with symbols and the dashed line) with our best-fitting line (solid line). In both panels, the shaded region shows the 
16th to 84th percentile scatter of the COCO+COLOR data. 


Table 2. Our best-fitting parameters for the C200 — M200 relation for the functional from of the Eqn. 


Co 

Cl 

C2 

C3 

c 4 

C5 

34.988 

-1.9841 

8.039 x 10" 2 

-1.777 x 10“ 3 

-1.4557 x 10“ 5 

7.34152 x 10“ 7 


To better quantify this discrepancy, we have fitted eqn. 03 to a 
combination of COCO+COLOR data. This ’COCO+COLOR compos¬ 
ite’ set was obtained by augmenting the COCO data for halo masses 
above 5 x 1O 11 /i _ 1 M 0 with the COLOR data and further supple¬ 
menting it at the low-mass end with the data obtained from the pro¬ 
file stacking. Our best fitting parameters are presented in the Tablepl 
and the corresponding fit is marked as a solid black line in Fig.pf 
To fix the asymptotic freedom of the fit at the low-mass tail, we 
have used the data points from Diemand, Moore & Stadel ( 2005[ l; 
|Ishiyama et al.| ( |2013| l and |Anderhalden & Diemand] < |20 1 3 [ ~ For 
haloes with 1 < M 200 /< 10 b our fitted relation pre¬ 
dicts concentrations that are ~ 10% lower than the SC 14 fit, while 
for even smaller haloes this tendency flips and our fit predicts con¬ 
centrations that are systematically higher. However, such a differ¬ 
ence has only small effects on the predicted boost factors for the 
radiation flux of DM annihilation. For completeness, we also com¬ 
pare the median COCO c(M 2 oo) relation with the model of Correa 
|et al.| ( |2015c| > (hereafter C15) which is based on the mass accre¬ 
tion history of haloes (see also |Correa et al.1 .2015a b; Ludlow et al. 
|2014[ >. The prediction of the C15 model for a WMAP7 cosmology 
is shown in the right-hand panel of Fig.[5]as the dashed-dotted line. 
This model agrees to better than 5% with our data for haloes more 
massive than ~ 1O 9 /i - 1 M0. At lower masses, the model under¬ 


predicts the concentration, which most likely reflects the fact that 
the C15 model was built for NFW profiles rather than for Einastio 
profiles. 

In Figure[6]we depict the time evolution of the concentration- 
mass relation. For 2 > 2, we find a flattening of the relation at the 
high mass end, with the flattening moving towards lower masses 
for higher redshifts. The same flattening is present also at 2 ~ 0, 
but we do not see it in the COCO data since the simulation does not 
resolve very high mass haloes. The concentration is related to the 
characteristic density of a halo, which in turn reflects the mean den¬ 
sity of the universe at the time when the central part of the halo has 
collapsed ( e.g. |Gao et al.|2008[|Ludlow et al.|2014j l. This naturally 
leads to the observed flattening of the c(M, 2 ) relation for very rare 
and massive objects, since these have assembled only recently and 
therefore share the same collapse time. We also find an evolution 
in the slope of the relation for lower mass haloes. This is in agree¬ 
ment with the well established picture according to which haloes 
are build-up hierarchically. Haloes continuously increase their mass 
with time, so the concentration at fixed halo mass is determined by 
different objects for each redshift bin. To better understand the time 
variation of the halo mass-concentration relation, we express the 
halo mass in terms of the mass variance, cr(M, 2 ) (see Eqn. |2|). 
The corresponding median c[cr^ 1 (M, 2 )] relation is given in Fig- 
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Figure 6. The redshift evolution of the concentration-mass relation for re¬ 
laxed haloes in the COCO volume. We plot the data from z = 9 (red colour) 
down to z = 0 (black). Solid lines mark the median Einasto concentration, 
while the filled regions show the lcr uncertainty in the relation. 



Figure 7. The halo concentration, C200, as a function of peak mass vari¬ 
ance, cr(M, z), for relaxed COCO haloes. We show the relation at various 
redshifts, from z = 9 (red) down to z = 0 (black). The black solid line 
gives the prediction of the Diemer & Kravtsov (2015 1 model at z = 0. 


ure [7] and shows that at fixed a(M) values the concentration varies 
only slowly with time. For comparison, we also give the prediction 
of the |Diemer & Kravtsov| |20 1 5} model to find that in the ‘big- 
peak’ regime (haloes corresponding to rare peaks at a given epoch) 
our results are in reasonable agreement with their prediction. At 
small <t -1 (M, z) values our data suggest a less steeper slope. The 
difference can be accounted for by noting that |Diemer & Kravtsov| 
( |2015j l have used a different halo finder and that their simulations 
had much lower mass resolution, hence they could not probe the 
’very small halo’ regime, which our simulation resolves. 



M fof [M 0 /h] 


Figure 8. The mean half-mass formation time, Z]/ 2 , as a function of present 
day halo mass, Mfof- The red points are the simulation data with error 
bars reflecting the bootstrap error on the mean. The shaded grey region 
marks the la scatter around the mean. The solid lines illustrate the fit to the 
COCO (red line) and to the MS-II (blue line) data. 


3.3 Formation times 


The highly hierarchical character of structure formation is espe¬ 
cially strongly imprinted in the build-up and mass accretion histo¬ 
ries of dark haloes. A significant fraction of a halo’s mass is as¬ 
sembled via the accretion of other haloes so it is natural to expect 
that more massive objects form later than the low mass ones, as 
confirmed by many N-body simulations of the CDM model. How¬ 
ever, the precise form of the halo mass-formation time relation and 
its intrinsic scatter is still a subject of discussionl |Lacey & Cole| 
|1993| [Wechsler et al.|2 002). Determining this relation is relevant 
for a number of reasons. Among which, most importantly, is un¬ 
derstanding how well this relation is correlated with the properties 
of the galaxies and the haloes they reside within, which is crucial 
for all models of galaxy formation (e.g. |Cole et al.|2000| Benson 


|et al.|2003]|Bower et al.|20061|De Lucia & Blaizot|2007||Guo et al. 

|20T3l >. 

The simplest and most used definition of the halo formation 
time, Zf, is given as the epoch at which a halo’s main progenitor 
assembles a fixed fraction of the present day halo mass. To compute 
this, we follow the main progenitor branch of each halo merger tree 
(see §2.3| for details on the merger trees) until the main progenitor 
reaches half of the halo’s final mass. Thus, our halo formation time 
is the redshift of half-mass assembly, Z 1 / 2 . This half-mass forma¬ 
tion time is one of the most commonly used formation time defini¬ 
tion found in the literature (however see |Li, Mo & Gao|2008| for 
other possible definitions). Figure [8] shows the formation time as 
a function of halo mass for our N-body simulations. For compari¬ 
son, we also show the fit to the relation obtained by|Boylan-Kolchin| 
|et al.| {2009 1 for MS-II haloes (rescaled from their M v i r to our mass 
definition Mfof ) and a fit to our own data. For the latter we use 
the same linear fitin log(l + Zi/ 2 ) asjBoylan-Kolchin et al. i2009 1 , 
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which takes the form: 


1 + Zl/2 — -<4-1 


= An ( - 




FOF 


\10 lo h- 1 MQ 


(13) 


The best fit to our data was obtained for Aq = 2.77 and /3 = 
—0.0765 and, as can be seen from the figure, it provides a very 
good description of the data over the entire mass range. 

Compared to the MS-II result we find that our haloes formed 
at similar epochs, however our formation time-mass relation is 
steeper. The difference can be presumably accounted for by the dif¬ 
ferent cosmology and due to different mass definitions. Most likely, 
the higher value of the erg used in the MS-II is a major driver of the 
discrepancy, as it is well established that this parameters is strongly 
correlated with the abundance of very massive haloes and hence 
also with the rate of their mass assembly. 

Figure [ 8 ] illustrates another important point. The halo-to-halo 
variation in formation times depends on halo mass, being the largest 
for low mass haloes and decreasing with increasing halo mass. This 
trend in halo-to-halo scatter can be understood in light of halo as¬ 
sembly bias, which was first pointed out by GaojfcJ¥hite|i 2007j(. 
The halo assembly bias (see also|Croton, Gao & White [2007 Lrj 
[Mo & Gao|2008) states that the halo formation time, Zf, is neg¬ 
atively correlated with the local amplitude of the density cluster¬ 
ing. It indicates that haloes form earlier in higher-density locations, 
which naturally have higher clustering amplitudes. Spatial regions 
characterised by a higher local amplitude of the density field lead 
to a faster DM halo formation for a number of reasons: (i) because 
of the excess matter clustering haloes can accrete more mass in the 
same unit of time, when compared to field haloes; (ii) higher clus¬ 
tering amplitudes imply also higher halo merger rates; (iii) the in¬ 
creased local density (compared to the universal background value) 
allows some density peaks to reach more rapidly the critical density 
threshold, S c , required for collapse. This reasoning can be inverted, 
when applied to regions with lower clustering amplitude than the 
mean, where exactly the opposite processes will induce later halo 
formation times. Another crucial ingredient is that DM haloes are 
biased tracers of the underlying DM density field (e.g. Frenk et ah] 


[T9851 |Davis et al.|i985] |Frenk et al.|l 9881 Col e & Kaiser|1989| >, 

with high-mass haloes having a large positive bias, while low mass 
haloes are slightly anti-biased (prefer to reside in lower-density en¬ 
vironments). Thus, massive haloes can only be found in regions 
characterised by a significant clustering excess. These regions are 


the nodes and the filaments of the cosmic web (see e.g. Bond, Kof- 


man & Pogosyan 1996 

Sheth & van de Weygaert|2004 Springel, 

Frenk & White 2006: 

Aragon-Calvo, van de Weygaert & Jones 


anisotropic nature of gravitational collapse. In contrast, the lower 
mass haloes pervade the whole range of large-scale environments, 
from voids to cosmic nodes, spanning many orders of magnitude 
in characteristic density (for the most recent results see e.g. Cautun 


let al.|2014c||Metuki et al.|2015||Falck et al.|2014[|Nuza et al.|2QI4| t 

A detailed investigation of the origin and properties of the halo as¬ 
sembly/cosmic web bias down to the smallest accessible halo mass 
is needed, in order to obtain a better physical understanding of this 
mechanism and its implication for halo properties and galaxy for¬ 
mation. In principle the COCO simulation set, owing to its resolu¬ 
tion and volume, is very-well suited for such studies. However, this 
venture is beyond the scope of this paper and we leave it for future 
work. 


4 DARK MATTER SUBHALOES 

In hierarchical CDM cosmologies, a significant fraction of halo 
mass growth takes place via the accretion of lower mass haloes, 
which results in a rich substructure of orbiting smaller DM clumps 
called subhaloes. The spatial distribution and abundance, kinematic 
and internal properties, and orbit parameters of these subhaloes are 
subject of intensive study in modern cosmology. Rendering a firm 
insight into the various physical properties of subhaloes plays a piv¬ 
otal role in linking the observed properties of Galactic satellites and 
dwarf galaxy population of the Local Group to the physical nature 
of dark matter. In this section we study the properties of the DM as 
a function of the mass of their host halo. 


4.1 Mass and velocity functions 


It is well known that due to discreteness of N-body simulations, 
effects like over-merging, two-body scattering, phase-space grain¬ 
ing and force softening will affect the internal properties of haloes 
and subhaloes that are close to the resolution limit of the simulation 
(see e.g .|Shandarin & Zeldovich|1989[ Klypin et al.|1999a||Power 


et al. 2003j jSpringelet al.|2008||Abel, Hahn & Kaehler|2012||Hahn 

Abel & Kaehler|2013) . Among others, these effects lower the maxi¬ 
mum velocity, Vmax, of small mass objects whose Rmax is compa¬ 
rable to the gravitational force softening of the simulation. We ap¬ 
ply the correction formula proposed by |Springel et al.| ( |2008| Eqn. 
( 10 ) therein), which, to a good approximation, under assumed per¬ 
fect circular orbits, accounts for this effect. While we have done 
so for all our haloes and subhaloes, we found that, due to our very 
high spatial resolution, the correction has negligible effects for the 
vast majority of objects with Vmax ^ lOkms -1 . The subhaloes 
below this Vmax limit are strongly affected by numerical effects, 
as we discuss in detail in Appendix [A] and are not considered in 
our analysis. 

Fig. El shows the mean cumulative number of subhaloes as a 
function of gL = M su t/Mn 00 , namely the gravitationally bound 
subhalo mass in the units of its parent host halo mass. The resolu¬ 
tion of our COCO run allows us to reliably identify substructure of 
various size in hosts down to a halo mass of ~ lO 9 /t _ 1 M 0 . To take 
full advantage of this, we bin host haloes according to their mass, 
M 200 , grouping them in five samples: 1 — 10 x 10 9 , 1 — 10 x 10 10 , 
1 - 10 x 10 11 , 1 - 10 x 10 12 and 1 - 1.5 x 1O 13 /i _ 1 M 0 . 

To keep consistency with previous works and in order to make 
fair comparisons, our analysis includes all subhaloes within a ra¬ 
dius, 7*50, from the host centre, as adopted in the analysis of the 
AQUARIUS suite (Springel et al .|2008| ). The radius, rso, is defined 
as the boundary at which the spherically averaged density reaches 
a value of 50 times the critical density for closure. On average, for 
galactic mass haloes, rso — 1-66 x r 2 oo, thus, one will find more 
subhaloes within r 50 than within r 2 oo- 

The AQUARIUS and Phoenix (Gao et al.j20l2| simulations 
have indicated that the substructure fractional mass function is well 
fitted over five orders of magnitude by a single power-law: 


N(> /r) oc fi s . 


(14) 


The best-fitting power-law suggests that s = 0.97 ± 0.02 for the 
Phoenix haloes and s = 0.94 ± 0.02 for the AQUARIUS suite. 
Both simulations have similar resolutions (for their highest level), 
but simulate host halo samples of different masses. The AQUARIUS 
hosts have an average mass of ~1O 12 /i _ 1 M 0 , while the Phoenix 
ones corresponds to a halo mass, M 200 a fewx 1 O 14 /i _ 1 M 0 . 

We fitted the same power-law to the substructure fractional 














































12 Wojciech A. Hellwing et al. 


Table 3. The best-fitting values for the power-law index, s (see Eqn.( |14| ), describing the subhalo mass function for hosts of different masses. 


(M200) in h *M 0 

10 10 

10 11 

10 12 

10 i: 

a 


10 13 

10 1 ' 

b 


s 

0.92 ± 0.02 

0.93 ±0.01 

0.94 ± 0.01 

0.94 ± 0.02 

0.95 ±0.01 

0.97 ±0.02 


a from AQUARIUS simulation 
^ from Phoenix simulation 
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Figure 9. Top panel: The mean cumulative number of subhaloes as a func¬ 
tion of p = M,. v i, /M 200 binned for different host masses. The coloured 
solid lines represent hosts from a halo mass of (M200) — 10 13 h~ 1 Mg 
(green) down {M 200 ) = 10®/i —1 Mq (brown). The black solid line shows 
the best-fitting power law function, N(> ji) ~ p,~ S ,withs = 0.95±0.01, 
found for the most massive host mass bin. The purple dashed line illustrates 
the scale-free mass function with s = 1. Bottom panel: The relative ratio 
of the data in each host mass bin to the COCO best-fitting power law model 
with s = 0.95. 


mass function for the COCO hosts, to obtain the best-fitting s pa¬ 
rameter as a function of host mass. The best fitting values and their 
standard error are given in Table [3] and were obtained by counting 
subhaloes with more than 100 particles. The best-fitting power-law 
function, for our best resolved hosts, which have a median mass, 
(M 200 ) = 10 13 /i _1 Mq, is shown as a solid line in Fig[9] Inter¬ 
estingly, this best-fitting value is found to be exactly in between 
the AQUARIUS and Phoenix results, with s = 0.95 ± 0.01, but 
consistent with those within the fit errors. Table [3] suggests that 
the power-law exponent, s, may increase very weakly with halo 
mass, but. given the error associated with s, this trend is not sta¬ 
tistically significant and a much larger study is needed to confirm 
or disprove such a trend. For a closer examination, we show in the 
lower panel of Fig. [9] the fractional difference with respect to the 


COCO best-fitting value of s = 0.95. The panel illustrates that for 
all mass-binned samples there is a range in fi for which the frac¬ 
tional difference exhibits an approximately flat region. At low /r, 
the deviations from a flat shape are driven by numerical resolution 
effects, while the behaviour observed at p, > 3 x 10~ 2 reflects the 
well known exponential cut-off in the mass function of the most 
massive substructures. 

The best-fitting power law exponents, s, that we found are 
close to the case of a scale-free subhalo mass function with the 
critical value of s = 1. For s = 1, each logarithmic bin in /r has an 
equal contribution to the total mass in subhaloes, which is logarith¬ 
mically divergent as /r —> 0. If the real substructure mass function 
is described by s = 1, than a significant fraction of the host mass is 
contained in subhaloes beyond the resolution limit of our simula¬ 
tion. For our best resolved sample with (M200} = 10 13 /i _1 Mq, an 
average of 7.7% of the host halo mass is contained in resolved sub¬ 
structure. Extrapolating this down to an Earth mass, corresponding 
to fi = 10 -13 , yields a fraction of 34% mass locked in substruc¬ 
ture. This prediction can be used further to yield DM annihilation 
gamma-ray flux (see e.g. |Bergstrom, Ullio & Buckley|1998}|Gon-| 
|dolo & Silk 1999). Detailed investigation of this situation is how¬ 
ever beyond the scope of this paper and we leave it for the future 
work. 

It has been long postulated that the scaled subhalo velocity 
function is independent of host halo mass, when expressed as a 
function of v = Vmax/V200 , which is the ratio between the sub¬ 
halo maximum circular velocity and the host virial velocity (e.g. 
|Moore et al.|1999| |Kravtsov et al.|2004) . Recently, this has been 
thoroughly confirmed using a large number of host haloes ( |Wang| 
|et al.|2012||Cautun et al.|2014b| >. Given both the very high resolu¬ 
tion and the large number of haloes in the COCO simulation, we can 
investigate the postulated invariance of the scaled subhalo velocity 
function over a wider dynamical range in subhalo V max and down 
to lower host halo masses. This is shown in Fig. [TO] where we plot 
the mean cumulative satellite count, N(> v), as a function of v for 
hosts binned according to their halo mass. To better highlight the 
invariance with host halo mass, the lower panel of Fig.|10|shows the 


|tun et a l. (2014b) for the mean subhalo count around galactic mass 
haloes. Since |Cautun et ak| ( |2014bf does not compute the subhalo 
count within rso, we take their result for subhaloes found within 
a distance of r 100 from the centre of the host halo. This leads to 
us counting more subhaloes than Cautun et al.]l |2014b| , which ex¬ 
plains why the COCO results are systematically above the zero level 
in the bottom panel of the figure. 

We find that the mean subhalo count, N(> v), exhibits at 
most a very weak dependence on host mass. This should be com¬ 
pared to the satellite abundances of Fig. [9] which show a strong 
dependence on host mass, M 200 • Any systematic deviations from a 
flat line for the results shown in the bottom panel of Fig. |10| appear 
only below the resolution limit of the simulation, which is show 


ratio between N(> u) measured in COCO and the best fit of 


Cau- 
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Figure 10. Upper panel: The scaled subhalo maximum circular velocity 
abundance as functions of the ratio v = Vmax/V 200 for hosts of different 
masses. The vertical arrows mark the v value corresponding to Vmax = 
10 kms —1 , which is our resolution limit. The black thin line shows the 
best-fitting function of |Cautun et al.|(2014b| for subhaloes within a radius, 
rioo, from the host. Lower panel: Ratios between the mean subhalo count 
at various host masses and the Cautun et al. result. 



V max [ km/s ] 

Figure 11. Upper panel: The subhalo count as a function of maximum cir¬ 
cular velocity for host haloes of different masses. The solid lines show the 
subhalo count using the present (z = 0) Vmax values. The dotted lines 
depict the subhalo count using the V max value at the subhalo's infall time. 
Lower panel: The ratio between the subhalo count using present time V max 
and the subhalo count using the Vmax at infall time. 


by a vertical arrow for each halo mass bin. The deviations seen at 
high v values for the most massive bin, A /200 = 10 13 /i _1 Mq, are 
due to the small number of host haloes present in that sample and 
hence are not significant. Thus, our results confirm the postulated 
invariance of N(> v) over four orders of magnitude in host mass, 
showing that this assumption holds for the majority of DM haloes 
that can host galaxies (IO 10 $ A/ 2 OOA/A /0 < 10 13 ). This invari¬ 
ance was exploited by Wang et al. (2012} and |Cautun et al.| i (20l4a| 
to derive new theoretical constrains on the mass of the Milky Way 
halo. 

In Fig.[TT]we show the mean subhalo count, N(> Vmax), as 
a function of subhalo maximum velocity, Vmax■ Describing sub¬ 
haloes in terms of the maximum circular velocity is more closely 
related to observations, since Vmax is more easily measured in ob¬ 
servations. The figure compares the subhalo abundance using the 
present day Vmax as well as the subhalo count as a function of 
the maximum circular velocity at subhalo’s infall time, Vmlx- The 
Vmax values are obtained by tracing the merger tree of each sub¬ 
halo and taking the peak value of Vmax throughout the history of 
the subhalo. For most practical applications, Vmlx is we ll approx¬ 
imated by the peak value of Vmax, since, once a halo falls into a 
more massive object, it becomes the subject of intensive tidal strip¬ 
ping and so its Vmax value is very likely to decrease rather than 
increase. Fig. m shows that at fixed subhalo size, i.e. fixed Vmax 
values, the abundance of objects roughly increase by an order of 
magnitude for each order of magnitude in host mass. This scal¬ 
ing is most pronounced for sufficiently small objects. This scal¬ 


ing breaks down for the most massive subhaloes, since the subhalo 
abundance changes its shape from a power-law to an exponential 
decline (see Fig. [9] and [TO}. Interestingly, a similar scaling is found 
also for the subhalo abundance as a function of VmJ x . This can 
readily be seen from the bottom panel Fig. [IT] that shows the ratio 
1Z = N(> Vmax)/N(> Vmax)- Both these functions are calcu¬ 
lated using the same objects, found at z = 0 inside a distance, rso, 
from their host, so their values are not influenced by the destruc¬ 
tion or accretion of new subhaloes. In other words, we expect that 
the ratio 1Z and its departure from unity are a good proxy for the 
efficiency of subhalo tidal striping at fixed subhalo circular veloc¬ 
ity. For small subhaloes with Vmax ~ 20kms _1 , we find that the 
ratio approaches 1Z ~ 0.35 for all host masses except for the low¬ 
est mass bin. For this lowest mass sample, (A/ 200 } = 10 9 A^ 1 Mq, 
COCO has enough resolution to identify only the most massive sub¬ 
haloes and it does not capture the power-law like regime of the 
subhalo abundance function. The convergence of the 1Z ratios to¬ 
wards a single values indicates that the efficiency of subhalo tidal 
stripping is comparable in hosts that differ by four orders of mag¬ 
nitude in mass, provided that the considered subhaloes are suffi¬ 
ciently small in comparison to their host. 


4.2 The radial distribution 


The radial distributions of subhaloes is also a subject of intensive 
study (e. g. |Gao et al.||2004| Diemand, Moore & Stadel||2004 De 
|Lucia et al.|2004[ Nagai & Kravtsov|2005||Wang, Frenk & Cooper 
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Figure 12. Radial number density profiles of subhaloes. Left panel: The radial distribution of subhaloes in galactic mass hosts, i.e (M200) — 10 12 /i _1 Mq. 
The lines with symbols corresponds to subhaloes of various masses, as expressed in terms p = M su b/M 200 , while the black solid line with crosses shows 
the results for all subhaloes. The bottom panel give the ratio of the various subhalo mass bins with respect to the distribution of all subhaloes. Right panel: 
The radial distribution of subhaloes for different host masses. The lines with symbols show the results of COCO. The dotted line (dashed-dotted) show the 
equivalent substructure profile from the AQUARIUS A-2 (A-3) run. The bottom panel shows the ratio of each data set to the reference case, which we take as 
(M200) = 10 13 /i 1 Mg . 


|M3| , since understanding how DM substructures are distributed 
inside their host haloes is important for several reasons. Among 
others, the radial distribution of subhaloes is instrumental in con¬ 
necting the observations of satellite galaxies with the properties 
of the background cosmology; it serves as input for many semi- 
analytical galaxy formation models; and it is important for strong 
lensing studies ( e.g. Mao & Schneider 1998^ Metcalf & Madau 


|200T] [Metcalf & Zhao||20021 |Kochanek & Da al||2004l |Xu et al" 

|20l5l >. 


Springel et al.jl|20Q8| have found that the radial distribution of 


subhaloes is independent of subhalo mass for at least five decades in 
mass (see Fig. 11 therein). We further investigate this finding, since 
the large number of host haloes of the COCO run allow for much 
better statistics. In addition, we further extended the analysis of 
Springel et al. by studying how the radial distribution of subhaloes 
varies with host mass. 

The left-hand panel of Fig. [T2] shows the dependence of the 
subhalo radial distribution on subhalo mass. The subhalo popula¬ 
tion is split according to their rescaled mass, fi = M su b/M 2 oo, 
following which, we stack the radial profiles of 53 host haloes 
whose median mass is (M200} = 10 12 /i _1 Mq. In the bottom- 
left panel we show the ratio of each subhalo mass sample with re¬ 
spect to the reference 'all subhaloes' line. We find a large degree 
of self-similarity between subhaloes of different masses, in agree¬ 
ment with the results ofiSpringel et al. |2008}. Flowever, we do 


find a weak, but systematic trend with subhalo mass. This is espe¬ 
cially pronounced for the most massive subhaloes with fj. > 10 -4 
for which, when compared to the distribution of all subhaloes, the 
radial distribution has an excess for r > r 2 oo and a scarcity at 
smaller radii. The exact size of this systematic effect is difficult 
to pinpoint because of the relatively large uncertainties associated 
with our data, which are caused by a significant host-to-host scatter 
in the distribution of massive subhaloes. 

The top-right panel gives the radial distribution of all sub¬ 
haloes for various host halo masses. In addition, we also show the 
results of the AQUARIUS A-2 and A-3 runs, to find that the COCO 
data for the best resolved haloes of mass (M200) = 10 13 h -1 Mg 
agrees with the AQUARIUS results down to a radial distance of 
~ 0.3r2oo- Below that radius, COCO contains fewer resolved sub¬ 
haloes than the higher resolution AQUARIUS runs. A similar be¬ 
haviour is seen when comparing the A-3 results to the A-2 ones, 
which is its higher resolution counterpart, and also when comparing 
the (M200) = 10 11 /i _1 Mg sample to the ( M200 ) = 10 13 h -1 Mg 
one. The systematic difference between COCO and AQUARIUS at 
r ~ O.6-R200 is likely a manifestation of the fact that the AQUAR¬ 
IUS sample contains a single halo and hence it is an indication of 
the object-to-object scatter. From the bottom-left panel of Fig. [12] 
which shows the ratio with respect to the reference (M200) = 
10 13 /i _1 Mg sample, we find that the radial distribution agrees 
remarkably well down to r ~ 0 . 2 r 2 oo for hosts spanning three 
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the values obtained by considering all resolved subhaloes in our 
simulation, while the solid blue curve indicates the same relation 
computed for field, isolated haloes. The mean Rmax value at fixed 
Vmax is a crude measure of central (sub)halo densities, as objects 
with the same Vmax, but larger (smaller) Rmax values are char¬ 
acterised by lower (higher) central densities. Fig. [T3] shows that at 
fixed Vmax, the mean Rmax values for haloes are 50% higher than 
for subhaloes. This can be easily seen as the blue solid line in the 
bottom panel of the figure. This behaviour reflects a well known 
fact, that subhaloes tend to have more concentrated density pro¬ 
files due to the effects of tidal stripping, which become significant 
once the subhaloes fell into their respective hosts (see e.g. Moore 


|et al.|1999}|Diemand, Kuhlen & Madau|2007|[Springel et al.|2008| > 

The tidal forces truncate a subhalo’s density profile by removing 
the mass that is only weakly gravitationally bound to the object. 
Since field haloes are rarely the subject of severe tidal forces, no 
such stripping takes place. An even more interesting find is the sys¬ 
tematic difference, at fixed Vmax, between the mean R ma x values 
characterizing subhaloes found in host haloes of different masses. 
Hence, satellites with the same Vmax values tend to have system¬ 
atically higher R ma x values when found in central haloes of lower 


5 CONCLUSIONS 


Figure 13. The Vmax — Rmax diagram for subhaloes and haloes. Upper 
panel: the relation between the maximum circular velocity and the radius 
at which it is attained for subhaloes found in hosts of various masses (solid 
lines with symbols and error bars). The solid black line shows this relation 
for all subhaloes in the simulation, while the solid black line illustrates this 
relation for field haloes. Lower panel: the ratios of the different samples 
with respect to the 'all subhaloes’ sample. 

decades in mass (10 11 — 10 13 /i _1 Mq). We can conclude that the 
spatial distribution of low mass satellites (i.e.with p <K 10 -3 ) has a 
universal shape across hosts of different masses. This is yet another 
example of the self-similar character of DM haloes in CDM cos¬ 
mologies. Detailed analysis of the subhalo radial density profiles is 
beyond the scope of our current paper and we leave it for a future 
work. 


4.3 The Vmax Rmax relation 

A fundamental structural property of subhaloes is the relation be¬ 
tween the maximum circular velocity, Vmax, and the radius cor¬ 
responding to this maximum, Rmax■ For DM dominated objects 
like spheroidal dwarf galaxies, the kinematics of the stellar com¬ 
ponent can be related to the underlying DM density profile via the 
Vmax — Rmax relation, which is the basis of numerous cosmolog¬ 
ical studies based on the stellar kinematics of Local Group dwarf 


galaxies (see e.g. 

Boylan-Kolchin. Bullock & Kaplinghat 

2011 

2012];|Sawala et al. 

2013.2016 

2015,2014; Wang et al. 2012 

Cau- 

tun et al. 2014a. Di Cintio et a 

. 2013; Zolotov et al. 2012; Arraki 

et al. 2014). While it is not our intention to have a detailed study 


of the Vmax — Rmax relation, we would like to add to the discus¬ 
sion by presenting an interesting finding. Namely, in Fig. [13] we 
plot the Vmax — Rmax diagram for both haloes and subhaloes. The 
lines with error bars mark the mean relation for subhaloes found 
in hosts of different masses. The solid black line corresponds to 


Since the establishment in the late 80s and 90s of the CDM model 
as the standard model for cosmic structure formation, it has become 
a subject of extensive tests and scrutiny. To understand the pro¬ 
cess of galaxy formation and evolution in a hierarchical ACDM 
cosmology, we need detailed knowledge of a multitude of physi¬ 
cal processes that act over an overwhelming range of scales. The 
formation and dynamical evolution of haloes and subhaloes, from 
tiny DM specks of Earth mass to the most massive gravitation¬ 
ally bound objects, together with highly non-linear and compli¬ 
cated baryonic processes set the framework in which galaxies live 
and evolve in our Universe. The constant development of observa¬ 
tional techniques is calling for an improvement in our theoretical 
modelling and understanding of the crucial phenomena involved. 
For this reason, we need simulations with ever growing resolution. 
However, we also need to model large enough cosmic volumes to 
obtain reliable statistics for various objects, from dwarf to giant 
galaxies. This is where simulations like the Copernicus Complexio 
play a pivotal role, since they have both a very high resolution and 
a large cosmological volume. In this paper we have introduced a 
new simulation, the COCO, that can reliable resolve all substructure 
down to a Vmax ~ 10 kms -1 in a cosmological relevant volume of 
~ 2.2 x 10 A h~ 3 Mpc 3 . This simulation is the first of its kind and 
is meant to be part of a whole series of intermediate zoom-in runs 
implementing both more ACDM cosmic volumes but also alter¬ 
native DM physics like Warm or Self-Interacting DM models (see 
also |Bose et al.|2016j l. 

The following is a summary of our main results: 

• The FOF mass function matches the ST and Reed formulas 
over seven orders of magnitude in halo mass at z = 0. However, 
for the intermediate redshift range of 2 < z < 0.5, the ST formula 
tends to over-predict the number of collapsed objects. 

• We have observed a departure of the c — M 200 relation from 
a single power-law at lower halo masses, in agreement with the 
results of |Sanchez-Conde & Prada| l |2014| . We give a best fit to 
the COCO data that reliably describes the concentration-mass re- 
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lation of relaxed haloes over six decades in halo mass 10 s ^ 
M 20 o/(/i" 1 M q ) ^ 10 14 . 

• We have probed the redshift evolution of the c— M 200 relation 
in the redshift interval, 9 ^ z ^ 0, to find that it is monotonic for 
small halo masses. 

• The hierarchical nature of halo formation processes is con¬ 
firmed for seven orders of magnitude in mass. The object-to-object 
scatter of the halo formation time depends on halo mass, with lower 
mass haloes showing a significantly larger scatter. This most likely 
is a manifestation of halo assembly bias, reflecting the multitude of 
environments in which low mass haloes are formed and evolve. 

• We have confirmed the power-law character of the subhalo 
mass function, N(> /.i) oc [a~ 3 , down to a rescaled subhalo 
mass, fj, = 10 -6 . For our best resolved hosts, with median halo 
mass, {M 200 ) = 10 13 ft -1 Mg, we find a power-law exponent, 
s = 0.95 ±0.01. 


Cautun 


• We find that the power-law exponent, s, depends on the host 
halo mass. It varies from s = 0.97 ± 0.02 for cluster mass haloes 
(Gao et al. 2012) to s = 0.92 ± 0.02 for lO lo /i - 1 M 0 haloes. 

• Our data confirms over a wider dynamical range in subhalo 
sizes and down to lower host masses that the mean subhalo abun¬ 
dance, N(> u), when expressed in terms of v = Vmax/V 200 , is 
to a very good approximation independent of host halo mass. The 
best-fitting results for N(> v), which were proposed by 
jet al.| l j2014bj , match our data down to our resolution limit. 

• The radial distribution of galactic subhaloes is nearly indepen¬ 
dent of subhalo mass, albeit with a very weak trend. Due to a large 
host-to-host scatter, this trend becomes visible only once we aver¬ 
age over a substantial number of host haloes. In addition, the radial 
distribution of subhaloes is nearly universal for hosts differing by 
three orders of magnitude in halo mass. 

• Finally, we have found that at fixed Vmax the mean Rmax 
values of subhaloes depend on the host halo mass, with lower mass 
hosts having subhaloes with higher Rmax values. This most likely 
reflects that at fixed subhalo size the tidal stripping processes are 
more efficient in more massive hosts. 


The current and future runs of the COCO suite will allow us 
to further test models of cosmic structure formation, including the 
development of semi-analytical galaxy formation models into the 
regime of low mass (sub)halo (hence also low galaxy luminosity) 
(for more details see |Guo et al.|[20r5| . This new satellite galaxy 
catalogue build on the base on COCO was alreay used for stringent 
statistical study of the prevalence of rare plannar sattellite configu¬ 
rations in the ACDM ( jCautun et al.|2 015). Moreover, our new set 
of simulations will allow for a better statistical study of radio-flux 
anomalies and lensing arc-distortions, the low-luminosity galaxy 
population, reionisation treatment in semi-analytical models and 
effects of the large-scale structures (Cosmic Web) on (sub)halo and 
galaxy properties and distributions. As these projects are currently 
work in progress, with the publication of this paper we also intend 
to make publicly available, in a short time, all relevant COCO halo 
and subhalo data bases accompanied by semi-analytical galaxy cat¬ 
alogues. In doing so, our hope and intention is to allow other re¬ 
searchers to use the COCO data for their own research projects. 
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APPENDIX A: NUMERICAL CONVERGENCE AND 
RESOLUTION TEST 


Here we assess a conservative limit for the mass and maximum 
circular velocity of the haloes and subhaloes that were resolved 
reliably in our simulations. This is necessary since haloes and es¬ 
pecially subhaloes that are resolved at low resolution are subject to 
many numerical artefacts that can alter their inner properties like 
density and circular velocity profiles, leading to unphysical results. 
In the following, we present two tests for determining the resolution 
limit of our numerical experiment. 

The first test investigates the relationship between the mass, 
M su b, and the maximum circular velocity, V ma x, of subhaloes. 
Following [Boylan-Kolchin et al.| |2010| l, we fit the relation with 
the power-law 


M sub = 6 x 10 1 


( V m ax 

\ 100 km s“ 


/i _1 M 0 . 


(Al) 


[Boylan-Kolchin et al.| have found that such a power law, with 
s = 3.23, provides a very good description of the MS-II data down 
to V m ax = 25 km s _1 . Below this value, the Vmax — M sub re¬ 
lation deviates from the power law fit suggesting that subhaloes 
with lower Vmax values are affected by numerical resolution ef¬ 
fects. The same power law, albeit with a slightly steeper scaling 
exponent, s = 3.3, gives a very good fit to the COCO data too, 
as shown in Fig. The COCO subhaloes follow the power law 
relation down to a much smaller V r 
5 x 10 7 /i -1 Mq). Fig. 


Al 


max values of ~ 10 km s 
also shows that that there is a 


(M a , 

very good convergence between the COLOR and the COCO runs, 
with COLOR having a resolution limit of Vmax = 25 km s“ . 
COLOR has the same behaviour as MS-II since both simulations 
have the same particle mass and force resolution. 

In the second test we compare the Vmax — Rmax relation of 
subhaloes with the same relation for haloes (see Fig. |A2j . We find 
that both for haloes and subhaloes, the median Vmax — Rmax re¬ 
lation shows an upturn indicative of numerical resolution effects at 
Vmax = 10 and 25 km s ^ 1 for COCO and COLOR, respectively. 
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Figure Al. A comparison of the Vmax — M su b relation for subhaloes 
found in hosts of all masses at redshift, z = 0. The COCO, COLOR and MS- 
II lines depict the median relation found by binning the subhaloes according 
to their Vmax values. The shaded green region illustrates the 16th to 84th 
percentiles around the median for the COCO sample. We also show the best 
fit power laws (Eq. EEU to the COCO and MS-II data. 



^max [km/s] 


Figure A2. A comparison of the Vmax — Rmax relation for haloes versus 
that for subhaloes. The solid yellow and orange lines show the relation for 
main haloes in the COCO and COLOR samples. The remaining lines show 
the Vmax — Rmax relation for subhaloes identified in all the COCO hosts 
(green solid line); in Milky-Way mass hosts (blue dashed line;see : |4.3} : and 
in the AQUARIUS level 2 hosts (solid brown line). The green shaded region 
gives the 16th to 84th percentiles for the full sample of COCO subhaloes. 


Thus, the two Vmax thresholds give a good conservative estimate 
for the resolution limit of our two simulations. 
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